First-principles study of the electronic and optical properties of Ho\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\text{W}}$$\end{document}W impurities in single-layer tungsten disulfide

The electronic and optical properties of single-layer (SL) tungsten disulfide (WS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2) in the presence of substitutional Holmium impurities (Ho\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\text{W}}$$\end{document}W) are studied. Although Ho is much larger than W, density functional theory (DFT) including spin-orbit coupling is used to show that Ho:SL WS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2 is stable. The magnetic moment of the Ho impurity is found to be 4.75\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _B$$\end{document}μB using spin-dependent DFT. The optical selection rules identified in the optical spectrum match exactly the optical selection rules derived by means of group theory. The presence of neutral Ho\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_W$$\end{document}W impurities gives rise to localized impurity states (LIS) with f-orbital character in the band structure. Using the Kubo-Greenwood formula and Kohn-Sham orbitals we obtain atom-like sharp transitions in the in-plane and out-of-plane components of the susceptibility tensor, Im\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _{\parallel }$$\end{document}χ‖ and Im\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _{\perp }$$\end{document}χ⊥. The optical resonances are in good agreement with experimental data.


Bandstructure
Numerical analysis. All Numerical calculations are carried out by using DFT and with the use of Perdew-Burke-Ernzerhof (PBE) generalized gradient (GGA) parametrization 19 for exchange-correlation functional. Fully relativistic noncollinear and spin polarized DFT calculations were performed as implemented in the Synopsis Atomistix Toolkit (ATK) 2021.06 20 . For Ho W impurity calculations, we consider a supercell consisting of eight unit cells along each crystalline-axis direction of the monolayer plane (i.e., 64 W atoms and 128 S atoms) and then replace a single W atom with Ho atom as shown in Fig. 1a). We consider a large supercell of edge length 25.22 Å, to fix the inter-impurity interactions. The point group of WS 2 with Ho W defect is D 3h . The periodic structure of the superlattice allows one to characterize the electron states by the bandstructure ε n (k) , where k is the vector in the first Brillouin zone of the superlattice and n enumerates different bands. The sampling of the Brillouin zone was done for a supercell with the equivalent of a 32×32× 1 Monkhorst-Pack k-point grid for the WS 2 primitive unit cell with a cutoff energy of 400 Ry. For all calculations, structure is first geometrically optimized with a force tolerance of 0.05 eV/Å. The formation energy for the Ho W impurity is calculated by means of the relation E tot [Ho W ] and E tot [host] are the total energy of the system with and without the impurity, respectively, n i is the number of added (n i > 0) or removed (n i < 0) species of atoms during the formation of the impurity. µ i 's are chemical potentials of the W and Ho atoms, which are estimated from their corresponding bulk forms. The small value of the formation energy E f [Ho W ] = 0.846 eV, which corresponds to an energy increase of E f [Ho W ]/64 = 13 meV with respect to the unit cell of WS 2 , indicates that the Ho W impurity in 8 × 8 × 1 WS 2 is thermodynamically stable. The thermodynamic stability can be shown by means of the cohesive energy per unit cell    Both cohesive energies are negative and therefore both pristine WS 2 and Ho:WS 2 are thermodynamically stable. In agreement with the above argument using the formation energy, pristine WS 2 is a little more stable than Ho:WS 2 . Therefore we can conclude that Ho:WS 2 at Ho doping concentrations of ≤ 1.56 % is thermodynamically stable. This is consistent with the experimental evidence that Er:MoS 2 is stable at an Er doping concentration of 3% 11 , and Er:WeSe 2 is stable at an Er/Yb co-doping concentration of 1.5% (1% Yb and 0.5% Er) 12 .
Given the difference in atomic radii of Ho (1.75 Å) and W (2.1 Å), relaxing the Ho:SL WS 2 structure and comparing it with the bulk structure of SL WS 2 indicates that local strain is introduced by the Ho W impurity in SL WS 2 . This is evident from the slight bond length distortion near the Ho impurity. The atomic configuration of Ho doped SL WS 2 preserves the D 3h symmetry with 2.609 Å and 3.29 Å being the Ho−S and Ho−W bond lengths, respectively. In pristine WS 2 we find bond lengths of 2.42 Å and 3.18 Å for W−S and W−W bond lengths respectively. Thus, we obtain a local distortion of 8.75%.
We first obtain the results for bandstructure and electric susceptibility for pristine WS 2 as shown in Fig. 1b,c, values of the band gap (1.64 eV) and splitting of the valence band edge (425 meV) due to SOC, are in good agreement with previously reported values [21][22][23][24] . The crystal structure of SL WS 2 is three atoms thin, where W atom is sandwiched in between two S atoms (S-W-S) via strong covalent bonds. Pristine SL WS 2 is invariant with respect to σ h reflection about the z = 0 (W) plane, where the z−axis is oriented perpendicular to the W plane of atoms. Therefore, electron states break down into two classes: even and odd, or symmetric and antisymmetric with respect to σ h . d-Orbitals of the W and p (t,b) -orbitals (t and b denoting the top and bottom layers) of the S atoms give the largest contribution to the conduction and valence band structure of SL WS 2 23,25 . Based on the σ h s y m m e t r y, t h e e v e n a n d o d d at om i c or bit a l s a re s p a n n e d by t h e b a s e s Using first principle studies [21][22][23] , it is known that valence and conduction bands are primarily made from d x 2 −y 2 , d xy and d z 2 of W atoms, which transform as E ′ 1 , E ′ 2 and A ′ irreducible representations (IRs) of the C 3h symmetry group at the K and K ′ points, in the absence of SOC ( Table 2). The presence of SOC couples the spin and orbital angular momenta, thereby requiring the consideration of the double-group IRs. Double-group IRs can be obtained by multiplying single-group IRs with E 1/2 as shown in Table 1, where E 1/2 is the 2D spin representation. The spin-orbit states for pristine SL WS 2 are shown in Fig. 6.
The isolated Ho atom has 11 4f electrons which are shielded by the outer 5s 2 5p 6 electrons. In addition, Ho is trivalent in nature i.e. when placed in a crystal, Ho has the tendency of losing 3 electrons (1 from the 4f and 2 from the 5s orbitals). To confirm the trivalent nature of Ho, we perform a Mullikan population (MP) analysis and calculate the atomic charges in different orbitals of Ho in SL WS 2 . The MP analysis shows that the 4f orbital of Ho is populated with an electron charge of 10.3e, whereas the 5s orbital is populated only with 0.372e, indicating a deficiency of 2.323e on Ho. Therefore to a sufficient approximation we consider the Ho impurity as Ho 3+ ions in SL WS 2 . The deficient electron charge of 2.323e on Ho is distributed in the entire supercell of Ho:SL WS 2 . For example, the 6 neighboring S atoms share an extra electron charge of 0.84e, while the 6 next-nearest neighbor W atoms share an extra electron charge of 0.072e, when compared with the pristine case.
The bandstructure of WS 2 with Ho W impurities is shown in Fig. 2. Regular electronic states within the valence or conduction bands are depicted by black lines while LIS ( f -orbitals of Ho) are depicted by blue lines. Some of Table 1. Character table of the group C 3h with ξ 3 = 1 . Two common notations are used for the IRs of the single and double group. The reduction of symmetry from D 3h to C 3h is accompanied by the  Fig. 7.
The Kramers theorem states that for every energy eigenstate of a time-reversal symmetric system with halfinteger total spin, there is at least one more eigenstate with the same energy. In other words, every energy level is at least doubly degenerate if it has half-integer spin. It can be seen that Kramers degeneracy, which is a consequence of time reversal symmetry, is broken for LIS in Ho W :SL WS 2 . In Ref. 26 it has been shown that presence of Ho Mo impurity leads to spin polarization and results in long range ferromagnetic coupling between local spins. The local magnetic moment of the Ho W impurity breaks the time reversal symmetry and lifts the Kramers degeneracy. In order to confirm that indeed the Ho W impurity in SL WS 2 contains a magnetic moment, DFT calculations are performed by using spin-polarized GGA method. The results are presented in Fig. 3, where we show that the Ho impurity has a magnetic moment of 4.75µ B . Our spin-polarized DFT calculations show that the exchange correlation potential leads to a spin splitting for Ho W :SL WS 2 . In Fig. 3b the isosurface plot for spin  www.nature.com/scientificreports/ density shows that the main contribution to the magnetism is due to the f -orbitals of Ho atom while the bulk states do not show any magnetic moment, in contrast to what has been observed in Ref. 26 , where a long-range magnetic interaction is seen. The reason is that we consider a much larger supercell of 8 × 8 × 1 , as opposed to their 4 × 4 × 1 supercell, resulting in a dilution of the impurity concentration that suppresses long-range magnetic interaction. The Ho W impurity preserves the σ h symmetry and therefore can be described by the group D 3h 27,28 , the irreducible representations (IRs) of which are shown in Table 2. The double-group IRs are obtained from the single-group IRs by taking the direct product with E 1/2 , as shown in Table 3. Figure 4 shows three examples of LIS that appear in the band structure. The 3-fold rotation symmetry C 3 of the impurity is clearly visible. Molecular orbital theory. A Ho W impurity inside WS 2 looks similar to an atom in an effective electrostatic ligand field created by its neighboring six sulphur atoms. In this approximation molecular orbital theory (MOT) can be used. To identify the LIS in the DOS, the projected density of states (PDOS) showing orbital contributions of individual atoms is shown in Fig. 2. In addition to the contribution from the f orbitals of Ho W , contributions from the p orbitals of the nearest neighboring S atoms and from the d orbitals of next-nearest neighboring W atoms are present. This means that in the Hilbert base spanned by

an LIS state can be represented by
where the real coefficients a i 's can be extracted from the PDOS shown in Fig. 2. Since admixture of orbitals is only allowed if they belong to the same IR, many coefficients are zero. The MOT diagram of pristine WS 2 can be found in Ref. 10 . The resulting eigenstates, identified by their IRs of D 3h , match the continuum states of the bands in WS 2 , as can be seen from Ref. 29 .  Table 3. Double-group IRs obtained from single-group IRs for D 3h group. www.nature.com/scientificreports/ Analyzing the PDOS, it becomes obvious that the Ho f orbitals couple to both the p orbitals of nearest neighboring S atoms and d orbitals of next-nearest neighboring W atoms. The resulting MOT diagram including Ho LIS is shown in Fig. 5. The orbital energy ordering can be determined by comparison with the PDOS shown in Fig. 2. The highest occupied molecular orbital (HOMO) is a E 1/2 spin-orbit state with an orbital A ′ 1 singlet state. The lowest unoccupied molecular orbital (LUMO) is a E 3/2 spin-orbit state with a E ′ orbital doublet state, which matches the PDOS in Fig. 2. Although the Ho atom with an average atomic radius of 1.75 Å is substantially larger than a W atom with an average atomic radius of 1.35 Å, DFT shows that the Ho W impurity is stable in the WS 2 host crystal. Because of the strong lattice distortions there are relatively strong hybridizations between the Ho f orbitals and the W d orbitals, as can be seen in the bandstructure in Fig. 2.

Optical response
Since the f-orbital contribution to the LIS is large, the optical spectrum exhibits narrow peaks, reminding of atom-like optical transitions. The relative dielectric functions ε r of various TMDs have been measured in Ref. 30 . We evaluate the matrix elements of the dielectric tensor in three dimensions ( i, j = x, y, z ) using the Kubo-Greenwood formula for the electric susceptibility where p j pq = �uk|p j |vk� is the dipole matrix element between Bloch states r |u k� = ψ uk (r) and r |v k� = ψ vk (r) , V the volume of the crystal, f the Fermi function, and Ŵ = 0.01 eV the broadening. A vacuum separation of a 3 = 20 Å has been chosen in order to suppress not only electron bonding but also electrostatic interactions. In this limit the Bloch functions are localized on SL WS 2 . Consequently, we can use the approximation (1/V ) k z → (1/�a 3 ) , where is the surface area of SL WS 2 . In this case χ = a 3 χ , which has the unit of length, is independent of the vacuum separation. Using this definition, we present the in plane χ and out of plane χ ⊥ components of the 3D susceptibility tensor for Ho W impurities in SL WS 2 in Fig. 7. We focus on transitions between states near the conduction and valence band edges and inside the band gap with resonance frequency ω uv = |ε u − ε v | , where ε u is the eigenenergy of the Bloch state ψ uk (r).
For pristine SL WS 2 the point group symmetry at the K and K' points is C 3h (see Table 1). A general result from group theory states that an optical transition is allowed by symmetry only if the direct product Ŵ(|vk�)⊗ Ŵ(p j ) ⊗ Ŵ(|uk�) contains Ŵ(I) in its decomposition in terms of a direct sum. Ŵ(I) denotes the IR for the identity, i.e., A ′ for C 3h . The in plane and out of plane components of p j vu must be considered individually because they transform according to different IRs of the point group. The resulting optical selection rules are shown in Fig. 6, which agree with the ones obtained from group theory shown in Table 4. These selection rules corroborate the difference between the in-plane and out-of-plane band gaps E g|| and E g⊥ , respectively, which can be seen in the in-plane and out-of-plane susceptibilities Im[χ � ](ω) and Im[χ ⊥ ](ω) , respectively [see Fig. 1]. We predicted this difference in Refs. 8,9 , which has later been experimentally confirmed 31,32 . This difference has also been verified theoretically by means of DFT calculations with GW correction and the solution of the Bethe-Salpeter equation for in-plane and out-of-plane excitons for similar crystals 33 . The band structure for pristine WS 2 in Fig. 1b shows in-plane and out-of-plane band gaps of E g|| = 1.64 eV and E g⊥ = 3.12 eV, respectively, which are in reasonable agreement with Refs. 34,35 . Alternatively, it is possible to use the conservation of angular momentum to derive the optical selection rules. In pristine SL WS 2 the C 3 rotational symmetry relaxes the atomic optical selection rules m = ±1 for σ ± transitions and m = 0 for π transitions to m = ±1 ± 3 for σ ± transitions and m = 0 ± 3 for π transitions, whereby an angular momentum mismatch of ±3 can be transferred to or from the crystal lattice. The resulting optical selection rules match the ones obtained above from group theory and are also shown in Fig. 6. When www.nature.com/scientificreports/ using the approximation of a two-band model described by a Dirac Hamiltonian for the conduction (CB) and valence (VB) bands, our selection rules match the ones shown in Ref. 18 . Given the point group symmetry of impurities in a crystal, the LIS transform according to its IRs. In the case of the Ho W impurity the point group symmetry is D 3h , its character table shown in Table 2. The identity for D 3h is A ′ 1 . Table 5 shows the the selection rules for electric dipole transitions for the IRs. Note that the electromagnetic field couples only to the orbital part of the Bloch states. Therefore the we need to consider only the orbital IRs of D 3h . Remarkably, we show in Table 6 that several optical transitions are in good agreement with available experimental data for optical transitions in Ho 3+ :YAG.

Conclusion
Our results of electronic and optical properties of Ho W impurities in SL WS 2 reveal LIS inside and near the band gap and atom-like sharp optical transitions both in χ and χ ⊥ . Therefore, we argue that REAs in TMDs are good candidates for spin qubits. Let us elaborate further.
Atom-like sharp optical transitions suggest that the decoherence time should be very long, which is one of the main criteria for a spin qubit. As mentioned in the introduction, by choosing a host material free of paramagnetic impurities and nuclear spins, it would be possible to substantially increase the spin coherence time of the impurity spin, in this case the spin of a Ho W impurity. Therefore, let us compare rare-earth atom spins in TMDs with other currently existing solid-state spin qubits: • W has a weak abundance of 14% of nuclear spin ( 183 W) and S has a negligibly small abundance of 0.8% of nuclear spin 3/2 ( 33 S). These can be removed by isotopic purification. In stark contrast to that, electron spin The term ±3 is due to the C 3 rotational symmetry of the lattice. These selection rules corroborate the difference between the in-plane and out-of-plane band gaps E g|| and E g⊥ , respectively. The Bloch states at the K and K' points transform according to the IRs of the C 3h point group.  Table 4. Electric dipole selection rules in C 3h symmetry. σ represents in plane transitions while π represents out of plane transitions. www.nature.com/scientificreports/ qubits in quantum dots made of GaAs suffer from hyperfine interaction. The issue is that Ga and As cannot be isotopically purified because all naturally abundant Ga and As isotopes have nuclear spins. Comparing to NV centers in diamond, N has a 99.6% abundance of nuclear spin 1 ( 14 N) and 0.4% of nuclear spin 1/2 ( 15 N). Therefore the nuclear spins of nitrogen cannot be removed by isotopic purification, either. In addition, P1 N impurities and surface spins are paramagnetic impurities that also lead to decoherence of an NV qubit. Consequently, we expect a much weaker decoherence of the Ho spin state. • The location in the direction perpendicular to the plane of the 2D material of the rare-earth impurities is accurate on the atomic level. By contrast, in 3D materials, such as GaAs and diamond, impurities and defects are spread throughout the 3D materials. Therefore we expect enhanced quantum sensing due to accurate distance to target atoms. • 2D materials have clean surfaces, in stark contrast to diamond that hosts dark P1 nitrogen impurities with nuclear spins and surface spins. The spin coherence time of shallow NV centers in diamond within 30 nm of the surface degrades drastically due to increased electric and magnetic noise. Diamond surfaces are difficult to be etched and polished in a controlled way due to diamond's hardness. 36 Ce 3+ in YAG exhibits electronic decoherence times of T 2 = 2 ms under dynamic decoupling 14 . This is relatively long considering that 27 Al, the only naturally occurring isotope, has a nuclear spin of 5/2. Liu et al. have recently performed the Deutsch-Jozsa quantum algorithm on the electron spin of a Ce 3+ ion in YAG by means of phase gates with an operation time of t op = 0.3 µs 37 . This would allow for N = T 2 /t op = 6.66 × 10 3 quantum operations.
In the case of Er 3+ impurities in CaWO 4 , the spin coherence time is T 2 = 23 ms, without isotopic purification 15 . In principle, this would allow for N = T 2 /t op = 6.66 × 10 4 quantum operations.
Since WS 2 can be isotopically purified to have zero nuclear spin, we expect even longer decoherence times and better performance with REAs in TMDs. In Ref. 38 a six-hour coherence time for optically addressable nuclear spins in Eu 3+ :Y 2 SiO 5 has been demonstrated experimentally. Based on this long coherence time a nuclear spin qubit in TMDs could allow for N > 7.2 × 10 10 quantum operations.
These advantages suggest that REAs embedded in 2D materials made of TMDs might be vastly superior to GaAs spin qubits and NV centers in diamond and could pave the way to realizing scalable quantum networks, scalable quantum computing, and ultrasensitive remote quantum sensing. Table 5. Electric Dipole selection rules in D 3h symmetry. σ represents in plane transitions while π represents out of plane transitions.